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ABSTRACT 

We propose a new multiscale method to calculate the amplitude of the gradient 
of the linear polarisation vector, |VP|, using a wavelet-based formalism. We demon¬ 
strate this method using a field of the Canadian Galactic Plane Survey (CGPS) and 
show that the filamentary structure typically seen in |VP| maps depends strongly on 
the instrumental resolution. Our analysis reveals that different networks of filaments 
are present on different angular scales. The wavelet formalism allows us to calculate 
the power spectrum of the fluctuations seen in |VP| and to determine the scaling 
behaviour of this quantity. The power spectrum is found to follow a power law with 
7 « 2.1. We identify a small drop in power between scales of 80 < Z < 300 arcmin, 
which corresponds well to the overlap in the u-v plane between the Effelsberg 100-m 
telescope and the DRAG 26-m telescope data. We suggest that this drop is due to 
undersampling present in the 26-m telescope data. In addition, the wavelet coefficient 
distributions show higher skewness on smaller scales than at larger scales. The spatial 
distribution of the outliers in the tails of these distributions creates a coherent subset 
of filaments correlated across multiple scales, which trace the sharpest changes in the 
polarisation vector P within the field. We suggest that these structures may be as¬ 
sociated with highly compressive shocks in the medium. The power spectrum of the 
field excluding these outliers shows a steeper power law with 7 « 2.5. 

Key words: ISM: general — ISM: structure — ISM: magnetic fields — radio con¬ 
tinuum: ISM — methods: statistical — techniques: image processing 


1 INTRODUCTION 


Previous power spectrum analysis of the amplitude of the 
polarisation intensity vector, |P| = \JQ'^ + where Q 
and U are the Stokes parameters, in the Galactic plane has 
shown evidence of large-scale structures in the Galactic mag¬ 
netic field ( [Haverkorn et aL]|2003| [Stutz et al.||2014) ). The 
power-law behaviour of these spectra is expected to be re¬ 
lated to the energy transfer from larger to smaller scales 
in the turbulent fluctuations of the magnetic field. Power- 
law variations as a function of Galactic latitude have also 



polarisation angle, 0 = (1/2) arctan(17/Q), which also show 
evidence of large-scale variations in the Galactic plane. Ac¬ 
cording to Landecker et al. (20101, large-scale variations of 


6 are probably associated with the large-scale features of 
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the magnetic field aligned with the spiral structure of the 
Galaxy. On the other hand, fluctuations of 0 on smaller 
scales could be explained by a turbulent Faraday screen in 
front of a uniform polarised background. However, such per¬ 
fect conditions are almost never satisfied in the interstellar 
medium (ISM) and most fluctuations seen in Stokes Q and 
U are probably due to a combination of Faraday rotation 
and intervening polarised emission along the line of sight 
in the Galactic plane. For these reasons, the interpretation 
of the direction and the amplitude of P when considered 
separately is very difficult. 

[Gaensler et al.| ( |2011[ ) proposed calculation of the am¬ 
plitude of the gradient of P, |VP|, as a new technique to 
measure variations of the vector P in the Q-U plane. It is 
defined as. 



This quantity can trace changes in both the direction and 
the amplitude of the vector P. Acting as an edge detec¬ 
tor in a map, the gradient of P highlights areas of sharp 
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change in the magnetic held and/or the free-electron den¬ 
sity, which are most likely due to turbulent huctuations or 
shock fronts in the ISM. Since |VP| is only sensitive to the 
smallest scales, it is not signihcantly affected by the loss of 
large-scale structure in interferometric data. On the other 
hand, one disadvantage of using the gradient is that it may 
enhance noise present in the data and the distribution of its 


amplitude depends on the telescope resolution (Burkhart 
et al.|2012| ). 


Variations of the emission probability distribution func¬ 
tion (PDF) width as a function of angnlar resolution and an- 
gnlar scale were measured in thermal dust emission, infrared 


emission as well as in |VP| structures (Froebrich & Rowles 
2010 Burkhart et al. 2012[ Robitaille et al.||2014 l. In gen¬ 


eral, this may be explained by small-scale high-extinction 
cores or structures that are not present on larger scales 
and do not create the large skewness typical of the lognor¬ 
mal distribution usually measured in star formation regions. 
The PDF of gas, dust column density and |VP| is also ex¬ 
pected to reflect the signature of physical processes occur¬ 
ring in the medium, e.g. turbulence, gravitational collapse 


or shocks (Ballesteros-Paredes et al. 2011 Burkhart et al. 


2012 Schneider et al. 20131. All those physical processes 


are scale dependent: shocks usually produce fine scale struc¬ 
tures, gravitational collapse depends on the local density of 
the gas and turbulence has the ability to transfer energy 
from large to smaller scales. Since the free-electron density 
and the magnetic field can also be affected by shocks and 
turbulence in the ISM, fluctuations in the polarisation inten¬ 
sity traced by |VP| should also be present on a broad range 


of scales. Previous studies using |VP| maps I^Gaensler et al. 


2011 Burkhart et al. 2012 lacobelli et al. 20141 concen¬ 


trated their analysis on small-scale fluctuations, primarily 
for two reasons: (1) the gradient of a map only samples the 
smallest scales and (2) many radio observations are made in- 
terferometrically and miss information on large-scale struc¬ 
tures since they do not completely sample the Fourier plane. 
In this paper, we propose a method to generalise such |VP| 
analysis to multiple scales using data where single-dish mea¬ 
surements are present. 

The calculation of the gradient as an edge detector in 
two-dimensional images has found multiple applications in 
different fields. Canny (19861 has shown, with image analy¬ 
sis methods for computer vision, that there is an uncertainty 
principle related to the detection and the localisation of a 
noisy step edge. In the presence of low signal to noise data, 
the precision of the localisation of edges must be traded by 
applying the gradient to a Gaussian-smoothed image. [Canny| 
(19861 also show that the hrst derivative of Gaussians of dif¬ 
ferent width can be used directly as a multiscale edge detec¬ 
tor. In a similar vein, Mallat & Hwang (1992 I generalised the 
method using wavelet transforms in a singularity detection 
algorithm applied to one and two-dimensional signals. Later, 
this generalised method, called the wavelet transform mod¬ 
ulus maxima (WTMM) was used by Arneodo et al. (20001 
as a multifractal analysis tool. 

In this work, we propose a similar technique based on 
wavelet analysis to calculate |VP| as a function of scale for 
data where single-dish measurements are present (a descrip¬ 
tion of these data is presented in Section]^. A description of 
the resolution effect on the calculation of |VP| is presented 
in Section the wavelet formalism is presented in Section 
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Figure 1. The amplitude of the polarisation intensity vector, |P|, 
calculated from Stokes Q and U maps of the CGPS data. 


1^ the formalism is tested on simulations in Section]^ appli¬ 
cation to real data and discussion are presented in Sections 
and Hand conclusions in Section 


2 OBSERVATIONS 


The analysis in this work is applied to a field of the Canadian 


Galactic Plane Survey (CGPS; Taylor et al. 20031 at 1420 
MHz on polarised data including Stokes parameters Q and 
U. For this survey, interferometric data were observed with 
the Synthesis Telescope at the Dominion Radio Astrophys- 


ical Observatory (DRAG; Landecker et al. 20001. All ob¬ 


servations were then completed with lower spatial frequen¬ 
cies from the Effelsberg 100-m telescope and the DRAG 
26-m Telescope ( Taylor et al.|[2003 Landecker et al.||2010 |. 
The chosen field is a combination of four mosaics from the 
survey available at the Canadian Astronomy Data Centre 
(CADcQ It covers ~ 8° in Galactic longitude and ~ 7° 
in Galactic latitude. The field is centred at I — 82.65° and 
b = 0.98°. The amplitude of the polarisation intensity vec¬ 
tor, |P|, is shown in Fig. CGPS data have an angular 
resolution of ~ 1 esc 5 aremin (where S is the declination) 
and a pixel size of 18 arcsec. 


3 RESOLUTION EFFECT ON |VP| 

One advantage of the calculation of |VP| is that the re¬ 
sults are not signihcantly affected by missing large-scale 
structures in interferometric data not completed with single¬ 
dish measurements. However, by sampling only the small¬ 
est scales, the gradient may enhance the noise in the data 
( Burkhart et al.||2012 |. Given the angular resolution of the 
CGPS maps and their pixel size, the synthesised beam is 
over-sampled by a factor of ~ 3.3 (Taylor et al.|20()3 1, which 


^ http://wwwl.cadc-ccda.hia-iha.nrc-cnrc.gc.ca 
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Figure 2. A one dimensional signal representing one row of pix¬ 
els located at I = —0.72° from the CGPS Stokes Q image. From 
the top to the bottom are shown the original signal with a super¬ 
posed smoothed version with a Gaussian filter having a standard 
deviation of 2® pixels (red line), the first derivative of the original 
signal and the first derivative of the smoothed version. 


means that, for these maps, the gradient is sensitive to vari¬ 
ations smaller than the synthesised beam. 

To visualise the resolution effect on the first derivative 
of a signal, the derivative of a one-dimensional function and 
its smoothed counterpart are shown in Fig. The signal 
represents one row of pixels at / = —0.72° from the CGPS 
Stokes Q image. To create a smoothed version of the signal, 
a Gaussian filter with a standard deviation of 2® pixels is 
convolved with the original signal. The smoothed counter¬ 
part is shown superposed on the original signal in the top 
panel of Fig. In the second panel, we show that in spite 
of obvious variations on larger scales in the one-dimensional 
signal, the first derivative of the original signal is much more 
sensitive to the smallest scale variation. On the other hand, 
the first derivative of the smoothed signal highlights vari¬ 
ations on larger scales which are independent of variations 
seen at smaller scale. 

Figure shows the spatial gradient of the linearly po¬ 
larized emission, |VP|, for the CGPS field. In the left panel, 
filamentary structures normally identified from the calcu¬ 
lation of |VP| can hardly be distinguished from variations 
at small scales associated with noise. However, the “honey¬ 
comb” noise variation pattern caused by the survey map¬ 
making becomes clearly visible. The right panel shows the 
gradient calculated from Q and U maps smoothed using a 
Gaussian beam with a standard deviation of 4 pixels. Since 
the first derivative is only sensitive to fluctuations larger 
than the synthesised beam, filamentary structures, similar 
to those initially presented by Gaensler et al. (20111, are 
now clearly visible. 

The effect of the map resolution on PDF moments of 
I VP I values was previously observed by [Burkhart et ah] 
120121. They used a Gaussian smoothing with FWHMs from 


3 to 9 pixels on Q and U maps and measured significant 
decreases in all of the first four moments of |VP| for simu¬ 
lations of turbulence with supersonic Mach numbers. Figure 
and [^ illustrate that smoothing maps of Q and U not 
only changes the distribution of |VP| values but can also 
significantly change its structures. In the next section, we 
propose a multiscale analysis technique in order to visualise 
and quantify changes in the distribution of |VP| maps as a 
function of scale. 


4 DOG WAVELET ANALYSIS 
4.1 Formalism 


The convolution of a Gaussian beam with maps of Q and 
U before the calculation of their gradient can give access to 
variations and sharp changes of |VP| at different angular 
resolution. By applying the technique on multiple angular 
scales, i.e. by gradually changing the Gaussian beam width 
convolved with maps Q and U, it is possible to extend the 
analysis of |VP| images in the spatial frequency domain. 


Following the work of Ganny (19861 and Mallat & 


Hwang (19921, wavelet transforms can be used as a basis for 


developing a multiscale edge detector analysis. The wavelet 
transform of a signal consists of the convolution of a set of 
functions, called daughter wavelets, each of which represents 
a scaled version of a mother wavelet. One class of wavelet 
functions called the Derivative of Gaussian (DoG) is defined 
by 


V>(®) = (-!)"“ (2) 

where 


4>{x) 



-kr 




( 3 ) 


The second order (m = 2) DoG wavelet represents the 
widely used “Mexican Hat” continuous wavelet. Even val¬ 
ues of m create symmetric functions which are appropriate 
for most general applications of wavelet transforms. Odd val¬ 
ues of m create asymmetric functions which are useful for 
revealing directional trends in data. They can also be used 
as edge detectors for structures present in an image. For the 
purpose of this analysis, the order of the mother wavelet will 
take the value of 1 or 3. 

In this section, polarised data are considered as two- 
dimensional functions f(x), where x is the vector position 
in a x-y plane. The continuous wavelet transform of f{x) 
with the DoG wavelet can be expressed as 


fihx) 


fl=l ^ J tplll ^ {x' — x)]f{x)(Ax' 
/2 = J- x)]f{x)d^x', 


( 4 ) 


where ijji{x,y) = d"^(f>{x,y)/dx’^ and xl) 2 {x,y) = 

d^(j>{x, y)ldy'^ (a three dimensional representation of these 
functions for m = 1 is shown in Fig. |^. All convolutions 
can be computed in the Fourier domain, which increases the 
speed of calculation. The function / represents the wavelet 
transform of / and I, the scaling factor of the wavelet. For 
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Figure 3. The spatial gradient of linearly polarized emission, |VP|, at the original resolution (left) and for the smoothed Stokes Q and 
U maps. Smoothed maps are produced with a convoluted Gaussian filter having a standard deviation of 2^ pixels. 



(a) (b) 

Figure 4. Analysing wavelets (a) —and (b) —'ip 2 . (The negative sign is for a better visualisation of the functions.) 


m = 1, it is interesting to note that equation is equiva¬ 
lent to the calculation of the gradient of f{x) smoothed by 
dilated versions of a Gaussian beam: 

hbx) (5) 

where G) is the convolution operation. From this point of 
view, the wavelet transform gives us a useful mathematical 
formalism on which a multiscaled version of |VP| can be 
defined. According to the statements above, |VP(l,a;)| can 
now be defined as: 


|VP(1, a:)| = ^m,x)\^ + \U{l,x)\\ (6) 

where, referring to equation Q, 


\Q{l,x)\ = -I- \Q2{1,x)\^, 

_ ^ ( 7 ) 

\U{l,x)\ = ^J\Uiil,x)\^ + \U 2 {l,x)\^. 


Since the work of Arneodo et al. 


the mathemat¬ 


ical formalism described by equations § to in has been 
usually associated with the WTMM method. In order to ex¬ 
tend the |VP| analysis to multiple angular scales, here we 
apply some components of this method to the CGPS polar¬ 
ization maps, in conjunction with a number of complemen¬ 
tary methods inspired by the wavelet analysis techniques, as 
the A-variance. However, a complete multifractal analysis, 
the original motivation for WTMM methods, is beyond the 
scope of this paper. 


4.2 Maxima chains 


Visually, the most interesting regions in maps of |VPj are 
areas showing the sharpest changes of the polarisation vector 
P. According to the WTMM method, one easy way to high¬ 
light these regions is by calculating the “maxima chains” of 
modulus values of the gradient. 

In addition to the magnitude of the polarisation gradi¬ 
ent, one can also calculate the argument or the direction of 
VP, at each position in a map iGaensler et al.|[2011l: 
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arg(VP) = 

( 8 ) 

where, 


Csign = sign 


dx dy ^ dx dy ) ' 


(9) 


Following this definition, modulus maxima are positions 
where |VP| is locally maximum in the direction of arg(VP). 
Thus, for every pixel on all scales, the argument of VP is 
calculated and the associated magnitude |VP| is compared 
with adjacent pixels having a similar arg(VP) value: ori¬ 
entations are divided into only six different directions to 
take into account the pixelisation effect. After an iterative 
process for the entire map, maxima positions should lie on 
connected “maxima chains”. Those chains allow us to vi¬ 
sualise locations where strong fluctuations in the electron 
density distribution and/or magnetic field strength occur. 
Chains are also useful to visualise coherent structures that 
are “connected” through multiple scales. 


4.3 Wavelet power sprectrum 


Similarly to the Fourier transform, the wavelet transform 
conserves the total energy of the original signal. This prop¬ 
erty can be defined following the generalisation of the 
Plancherel identity for the continuous wavelet transform 
| Farge|1992 Dumas et al.||20fT l: 


J \fix)\^d'^x = J j 




( 10 ) 


where 


= 


/ 




( 11 ) 


Equation ( |11[ ) is also called the admissibility condition of the 
wavelet, where iplk) is the Fourier transform of the mother 
wavelet tp{x) and k is the wavenumber vector. This condi¬ 
tion is satisfied for every m-th order of equation H- From 
equation (101, the energy conservation can also be defined 


as a function of spatial scale only: 


Em- 1 




( 12 ) 


This relation shows that wavelet coefficients can be com¬ 
pared to Fourier coefficients and, as for the calculation of 
the Fourier power spectrum, wavelet coefficients can be used 
to measure the energy transfer from large to smaller scales. 
It is important to note that the normalisation factor 
in equation ( [^ is only required to ensure the validity of 
equation ||^ ( [Arneodo et al.|2000 1. In order to calculate the 
wavelet power spectrum of | VP(Z, x)\, the regular normalisa¬ 
tion of a wavelet transform, , is used. The wavelet energy 


spectrum defined by equation (121 can also be expres sed in 
terms of the Fourier energy spectrum, E{k) — l/(fc)|^ (Farge 


1992): 


E{1) = J E{k)\i>{lk)\^(fk. (13) 

This relation means that at a particular scale, the global 
wavelet energy corresponds to the integral of the Fourier 
energy spectrum of the analysed function weighted with the 
energy spectrum of the wavelet at that scale. 

In order to produce a wavelet power spectrum similar 
to the classical Fourier power spectrum, which takes into 
account the finite size of map f{x) and the discrete number 
of pixels, we use the relation: 

= (14) 

^ X 

The notation Sp{l) is used here instead of the usual P{1) 
for the power spectrum, in order to avoid possible confusion 
with the polarisation intensity P. 


4.4 Equivalence with Fourier wavelength 


If one wants to compare the wavelet scale I with the 
wavenumber k in the Fourier domain, the equivalence of the 
scaling factor I in the frequency domain has to be defined. 
The wavelet analysis described in the previous sections is 
equivalent to the calculation of the gradient of an image 
smoothed by Gaussian filters of different widths (see equa¬ 
tion §). Following this statement, the scale I defined in the 
previous sections is related to the standard deviation of the 
Gaussian filter. For the following analysis, the wavelet scal¬ 
ing factor will be defined as Zp = Z • (27r)“^, so that the scale 
Z can be directly compared to the wavenumber k = 1/1 in 
the Fourier domain. 

The mth derivative of the Gaussian makes the function 
oscillate around zero. The wavelength of these oscillations 
and their amplitude, rapidly decaying towards infinity, are 
the two properties which allow the wavelet function to be 
localised in the frequency domain. The width of the function 
in the frequency domain acts as a bandpass filter. According 
to Kirby (20051, one easy way to define the relationship 


between the scale of the wavelet function and the frequency 
content of the signal is to determine the wavenumber at 
which the wavelet function is maximum in the frequency 
domain. He determined that, in the case of the DoG wavelet, 
the equivalence between the wavelet scale and the wavelet 
scaling factor is Zs = y^mj/Zp, so that the scaling factor 
becomes Zp = ^/m ■ I ■ . 

In other words, for this analysis, the wavelet scaling fac¬ 
tor is not chosen to correspond to the standard deviation of 
the initial Gaussian, but it is chosen to correspond instead 
to the Fourier wavenumber k that is sampled by the band¬ 
pass filter in the Fourier domain. This definition allows a 
better comparison between the wavelet power spectrum and 
the classical Fourier power spectrum. 


5 TESTS ON SIMULATED DATA 

The formalism presented in the previous section shares sim¬ 
ilarities with the A-variance introduced by |Stutzki et al.| 
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(19981. This method has been successfully applied in sev¬ 


eral studies in order to characterise structures at multiple 
scales induced by turbulence in molecular clouds ([Bensch 


et al.| 200 T||Qssenkopf et al.|2008a|b| [Alves de Oliveira et ah 


2014). The A-variance is defined as a measure of the amount 


of structure at a given scale / in a map. Its definition is sim¬ 


ilar to the energy spectrum defined in equations (121 and 


(131, except that the convolved filter 'ip{x) is isotropic. As 
mentioned by Ossenkopf et al. ( 2008a| ), the main advantage 
of the A-variance method comes from its smooth filter shape 
which ensures a robust angular average of the signal and a 
lower sensitivity to singular variations and finite map size 


effects. Similarly to the work of Bensch et al. (20011 which 


tested the influence of telescope beam and hnite map sizes 
on the A-variance, this section tests those effects using the 
anisotropic DoG wavelet. 

The robustness of the wavelet power spectrum calcula¬ 
tion is tested on two Gaussian random field (Grf) simula¬ 
tions of 1024 X 1024 pixels for both Stokes Q and U maps. 
Those images are produced by applying a power law as a 
function of scale to the squared amplitude of a random- 
phase map. Similarities between Grfs and interstellar struc¬ 


tures was pointed out by Stutzki et al. (1998 1 in their study 


of fractal properties of molecular clouds. A Grf simulation 
with a power law of 7 = 2.5 representing Stokes Q maps is 
displayed in Fig.[^ The original Grf is displayed in Fig.j^a), 
Fig. §b) shows the same field convolved with a Gaussian 
filter having a standard deviation of 2 pixels and Fig. |^c) 
shows the original field added with a random noise having a 
Crms of 0.5. The original field has a mean pixel value of zero 
and a standard deviation of 1.0. The Fourier power spectra 
of |P| for the three fields are shown in Fig. (a). They are 
calculated on 2048 x 2048 extended maps with zero-padding 
and an apodised interface between the extension and the 
image on 5 per cent of the border of the original image. 
To avoid spurious power at smaller scales caused by edges 
of the image, the mean pixel value of images must be sub¬ 
tracted before the apodisation. In Fig.[^(a), the spectra are 
produced by averaging the squared amplitude of complex 
Fourier coefficients over different annuli in the u-v plane. 
Figure]^ (a) shows that, at small scales, the telescope beam 
and noise induce a significant departure from the power law. 

The wavelet power spectra of |VP| calculated follow¬ 
ing equations <§ and ( |14[ ) are displayed in Fig. (b) and 
(c). The calculation is done on the same extended maps 
as for the Fourier power spectra. Equidistant values of I in 
logarithmic scale are chosen starting from 4 to 1024 pixels 
with an interval of 2°'^® pixels. The same pixel resolution 
as the GGPS data has been assigned to the Grfs. Wavelet 
power spectra on figures [^(b) and (c) are respectively asso¬ 
ciated with the first and the third derivative of a Gaussian, 
i.e. m = 1 and 3 in equation ([^. A s noticed by Arneodo 
|et al.| ( |2000[ ) and |Khalil et ah] ( |2006| ) , the robustness of the 
results can be tested using the DoG wavelet by repeating 
the analysis with a wavelet of a higher order. The first and 
the third derivative of a Gaussian have respectively one and 
three vanishing moments. A wavelet with more vanishing 


^ The apodisation consists of the multiplication of a taper func¬ 
tion, in this case the negative slope of a cosine, which smooths 


moments can represent more complex functions. For exam¬ 
ple, the wavelet transform of a polynomial function of degree 
n will be equal to zero, if a wavelet has vanishing moments 
up to the order m ^ n ( Sheng|2010 1 . Gonsequently, repeat¬ 
ing the analysis with a wavelet of a higher order can confirm 
that the scaling behaviour of the wavelet transform is not 
dominated by the order of the analysing wavelet and can 
also highlight the effect of polynomial distributions changing 
the self-similar geometry of the data. For both orders of the 
DoG wavelet, a power law of 7 = 2.5 is fitted for 4 < Z < 50 
arcmin on the wavelet power spectra of the original Grfs. 
The wavelet with more vanishing moments is significantly 
more sensitive to the beam smoothing effect. The third or¬ 
der wavelet is also more affected by the noise, but less than 
by the beam convolution. It is important to note that the 
noise wavelet power spectrum with the hrst and the third 
order DoG wavelet has a flat power law, i.e. 7 = 0 instead 
of 7 = 2 as with the A-variance. This difference comes from 
the normalisation choice discussed previously in section [T^ 
The third order wavelet is also less affected by the edge ef¬ 
fect at larger scales than the first order wavelet. As shown 
in Fig. since wavelet functions decay as a:“", where n is 
the order of the wavelet, the third order wavelet has a bet¬ 
ter localisation in the spatial domain than the first order. 
For this reason, the third order wavelet is less affected by 
the zero-padding which decreases the power of large-scale 
structures. 

Figure (d) shows the values of the fitted power laws 
to the wavelet power spectra between 4 < Z < 50 arcmin, 
for five different power law indices of the original Grfs. For 
the first order wavelet analysis, an underestimation of the 
spectral index is measured for 7 > 3. This effect was also no¬ 


ticed by Bensch et al. (2001 1 for the A-variance analysis and 
was attributed to the fact that edge effects are significant for 
maps covering only a fraction of the spatially extended emis¬ 
sion. This statement is true for steep spectral index, where 
large-scale structures dominate. However, an overestimation 
of the spectral index is measured for 7^3 using the third 
order wavelet, even if edge effects are less important for this 
wavelet. In that case the overestimation of the spectral in¬ 
dex can be attributed to the lower resolution of this wavelet 
in the spectral domain. 

As for the A-variance, our wavelet power spectrum can 
satisfactorily recover the power law index of the fractal sim¬ 
ulations for 2.0 ^ 7 ^ 3.0. According to their statistical 
properties, Grfs have the same power law index in every di¬ 
rection. Following this property, the calculation of the power 
spectrum of |VP| using equations 1^ [tI and 


14 


the sharp edges in an image. 


can recover 

the power law index of individual Stokes Q and U simu¬ 
lated maps. Because of the normalisation choice discussed 
previously in section |4.3| the slope of the wavelet power 
spectrum of |VP| is equal to the power law of the Grfs. 
This is similar to the A-variance where the slope a is re¬ 
lated to the power law 7 following the relation a = I 7 I — 2 . 
Real Stokes Q and U data are spatially correlated but are 
not assumed to have exactly the same power law index. In¬ 
tervening polarised emission and faraday rotation along the 
line-of-sight should induce spatial correlation between Q and 
U maps and should also modify the measured power law of 
the wavelet power spectrum of | VP(Z, a;)|. Gonsequently, the 
wavelet power spectrum of |VP(Z,a;)| is a unique measure 
of the variations of the polarisation vector P in the Q-U 
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Figure 5. The Grf simulation of 1024 x 1024 pixels of Stokes Q map with 7 = 2.5: (a) the original Grf with a mean pixel value of zero 
and a standard deviation of 1.0, (b) the same field convolved with a Gaussian filter with a standard deviation of 2 pixels and (c) the 
original field with random noise, with an rms of 0.5, added. 
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Figure 6. Power spectra analysis of Grfs simulated Stokes Q and U maps: (a) The Fourier power spectrum of |P| of the original image 
(diamond), the image convolved with a Gaussian beam (star) and the image with added noise (plus), (b) The wavelet power spectrum 
of the three same images using the first order wavelet, (c) The wavelet power spectrum of the three same images using the third order 
wavelet, (d) shows the values of the fitted power laws to the wavelet power spectra (for 4 < / < 50 arcmin), for five different power law 
indices of the original Grfs. Diamonds represent power laws measured with the first order Dog wavelet (m = 1) and triangles represent 
power laws measured with the third order Dog wavelet (m = 3). 
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Figure 7. The positive part of the first and the third derivative of 
a Gaussian centred at zero, respectively m = 1 and 3 in equation 


§. 


plane as a function of the angular scale and should not be 
directly compared with the Fourier power spectrum of |P| 
or of Stokes Q and U maps alone. 


6 APPLICATION TO CGPS DATA 


The wavelet analysis technique described in Section was 
applied to the CGPS field shown in Fig. [^for both Stokes 
parameters Q and U. Each map was extended to the closest 
power of 2 - in this case, 2048 x 2048, with zero-padding 
pixels and apodised on 5 per cent of the border of the original 
image. Angular scales I are chosen following the same rules 
as for the Grf simulations. 

The amplitude of the gradient of P for the CGPS field 
at four different scales, overplotted with maxima chains, is 
displayed in Fig.[^ For clarity, maxima chains of the small¬ 
est scale, I — 9.6 arcmin, are only displayed for pixel values 
above 0.15 K (« Surms). Each wavelet transform in Fig. 

shows very different filamentary structures. The complex 
network of filaments at smaller scales is replaced by a more 
extended network of filaments on larger scales. The general 
pattern of the lower scale is sometimes preserved and some¬ 
times not. Particularly, some features described as “double 


jump” profiles by Burkhart et al. (20121 (see green boxes 
in left and right upper panels of Fig. appear only for a 
small range of scales. Such features are associated with the 
derivative of a delta function that can result from interac¬ 
tions of strong shocks. On the other hand, other subsets of 
the network persist over multiple scales and create a subset 
of “coherent” structures across the field. An example of a 
“coherent” subset network is displayed in Fig. The sub¬ 
set is selected using an iterative algorithm called the scale- 
wise Goherent Vorticity Extraction (GVE) (see [Nguyen van 
yen et ah] (|2012 \ and Robitaille et al. (20141 for a detailed 


description). As a function of scale I, this algorithm con¬ 
verges to an optimal threshold value to separate outliers, 


i.e. non-Gaussianitie^ from randomly distributed wavelet 
coefficients of |VP(Z,a:)|. Figurej^shows maxima chains for 
which the maximum value along the chain is part of the 
separated outliers. 

The wavelet power spectrum calculated according to 
equation (14 1 is shown in black diamonds in Fig. 10 and 11 


Wavelet coefficients |VP(Z, x) \ are calculated using the third 
order wavelet in order to highlight effects produced by noise 
and edges on the map (see section . The flattening and 
the power drop caused by the Gaussian apodisation func¬ 


tion applied to the DRAG Synthesis Telescope data (Lan- 
decker et al.||2010 Stutz et al.||2014 l clearly appears in the 
power spectrum between 1 < ^ < 6 arcmin. A flattening of 
the spectrum is also produced by the finite size of the map 
around 600 arcmin. The power spectrum shows a power law 
behaviour between 10 < / < 80 arcmin followed by a small 
drop of power between 80 < Z < 300 arcmin. This drop cor¬ 
responds well to the overlap in the u-v plane between data 
from the Effelsberg 100-m telescope and the DRAG 26-m 


telescope, which is between baselines of 3 to 15 m (Lan- 
decker et al.||2010 |. The corresponding angular sizes, 48 to 
240 arcmin, are indicated by dotted lines in Fig. |10[ The 
26-m data were initially undersampled and gaps where no 
26-m data were available were filled with smoothed Effels- 
berg data. This undersampling, and the process applied to 
correct the data, could have produced an underestimation of 
the power over that range of scales. The average calibration 
ratio between the two datasets is 0.96±0.01 (26-m/lOO-m). 
This small factor could explain the drop in power seen in the 
wavelet power spectrum of the GGPS field. For this reason, 
a power law is fitted only between 10 to 60 arcmin following 
the relation Sp{l) = So ■ P, where the fitted values of pa¬ 
rameters So and 7 are (4.69 ± 0.03) x 10“^ and 2.15 ± 0.01 
respectively. 

The wavelet power spectrum has also been calculated 
using only the Gaussian coefficients of \\7P{l,x)\ separated 
by the scale-wise GVE. This Gaussian power spectrum is 
represented by the red stars in Fig. [m Both distributions, 
for all coefficients and for the separated Gaussian part, are 
plotted with the black lines in Fig. |12[ The distributions are 
normalised following the definition: 


I{l,x) = 


(|VP(/,a:)|), 


(15) 


where {}x, is the average operator over all x. The normalised 
distributions for all coefficients between scales of 6.8 to 91.3 


arcmin on the top panel of Fig. 12 (black lines) are lognor¬ 


mal and consequently, the average value of the coefficients 
does not accurately characterise the distribution. The sep¬ 
arated Gaussian part shows a peak centred on I{l,x) « 1, 
which means that the average value of coefficients is more 
representative of the general tendency of the distribution. 
The fitted parameters for the Gaussian power spectrum are 
(6.4 ± 0.5) X 10“® and 2.52 ± 0.02 for So and 7 respectively 
for scales between 20 to 60 arcmin. Scales 108.6 < / < 614.4 


^ By construction, following eq.the distribution of P{l,x)\ 
cannot respect a perfect Gaussian distribution, however the ter¬ 
minology Gaussian and non-Gaussian are used to describe respec¬ 
tively the symmetrical part and the tail of the distribution. 
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Figure 8. From left to right: in “cubehelix” colour scale I Green|[2011 1 are the |VP(Z,a;)| values for the same field as in Fig. at four 
different scales 1= 9.6, 45.7, 153.6 and 434.4 arcmin. White lines represents maxima chains (see section [4.2[| corresponding to the scale. 


arcmin do not respect the self-similarity of small-scale distri¬ 
butions. However, a clear separation between two different 
behaviours as a function of scale is hard to establish and the 
transition from lognormal to Gaussian distributions might 
also be continuous. 


7 DISCUSSION 

The new method proposed in this paper allows us to extend 
to multiple scales the study of structures produced by the 
calculation of |VP|. As shown in Fig. filamentary struc¬ 
tures in I VP(/, x)\ are highly dependant on the angular scale 
or on the instrumental resolution with which the polarised 
signal is observed. Furthermore, the wavelet power spectrum 
analysis of \\7P{l,x)\ allows us to identify scales where the 
signal to noise ratio or the beam signature becomes impor¬ 
tant and causes a flattening of the power law behaviour. 
Consequently, studies of the gradient of P applied only to 
the smallest spacial scales, without or with little smoothing 
of original data, should be aware that a significant amount 


of the structure seen at lower intensity in |VP| images may 
be associated with noise. 

Turbulence is expected to be one of the major pro¬ 
cesses responsible for fluctuations on multiple scales in the 
ISM. The power law measured across a large range of scales 
(10 ~ ^ ~ 60 arcmin), could be associated with the pres¬ 
ence of turbulence in the magnetic field. The self-similarity 
of wavelet coefficient distributions plotted in Fig. is an¬ 
other indication that turbulence can play a major role in 
fluctuations seen in the magnetic field and/or the electron 
density over the same range of scales. Although, as demon¬ 
strated with magnetohydrodynamic (MHD) simulations by 
|Gaensler et al.| | |2011[ ) and [Burkhart et al.| ( |2012| ) , various 
types of turbulent environment, e.g. subsonic or supersonic 
turbulence with different Mach numbers, can create various 
type of structures in |VP| maps. As expected for a large 
field localised in the Galactic plane, the filament network of 
I VP I displays a large range of intensities and different types 
of discontinuity associated with different types of fluctuation 
in the magnetic field and/or the electron density. We see in 
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Figure 9. The superposition of maxima chains from scale I =22.8 to 258.3 arcmin over the map of IVPI at Z = 22.8 arcmin. They 
represent a subset maxima chains for which the maximum value along the chain is part of outliers separated with the scale-wise CVE 
algorithm. 



Figure 10. The wavelet power spectrum of |VP(Z,cc)| for the 
CGPS field for all coefficients (black diamonds). 


Fig. that even if the wavelet coefficient distributions are 
self-similar for a given range of scales, the power is not ran¬ 
domly distributed. The tails of the lognormal distributions 
shown in the upper panel are associated with coefficients 
that contribute more to the average power at a given scale 
than randomly distributed coefficients. This excess has a 
significant impact on the measured power law. The Gaus¬ 
sian part of the distributions for the self-similar range of 
scales possess a steeper power law (7 « 2.5) than the orig¬ 
inal distribution (7 « 2.1). The spatial distribution of non- 
Gaussianities displayed for a small range of scales in Fig. 
also shows that these structures tend to create a coherent 
subset of filaments correlated across several scales. These co¬ 
herent structures can have different origins. [Burkhart et al.| 
| |2012| ) showed that moments of |VP|, i.e. mean, variance, 
skewness and kurtosis, are correlated with the Mach num¬ 


Figure 11. The wavelet power spectrum of |VP(Z,cc)| for the 
CGPS field for all coefficients (black diamonds) and for the Gaus¬ 
sian part of the distribution (red stars). The solid line represents 
the Fourier power spectrum of |P|. 


ber of MHD simulations. Higher Mach numbers create more 
asymmetric distributions which have tails at high intensity. 
Higher intensity structures in |VP| associated with those 
tails are caused by sharp changes of the polarisation vector 
P that can be produced by compressive shocks in a super¬ 
sonic turbulent medium. In dense regions, the magnetic field 
lines are frozen into the ionised gas and compressive shocks 
will amplify the magnetic field intensity. Under these condi¬ 
tions, the magnetic field intensity is correlated with the elec¬ 
tron density and creates higher intensity structures seen in 
I VP I ( Burkhart et al.||2012 b However, subsonic turbulence 
induces no compressive motion and in that case fluctuations 
traced by |VP| are dominated by random fluctuations in 
the gradient of the magnetic field. 
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Figure 12. Normalised distributions of wavelet coefficients of 
|V-P{/,£c)| for the CGPS field for all coefficients (top panel). The 
black lines represent scales between 6.8 and 91.3 arcmin and the 
blue lines present scales between 108.6 and 614.4 arcmin. The 
lower panel shows the Gaussian part of the distribution for scales 
between 6.8 and 91.3 arcmin. 


It is interesting to note in Fig. that the power spec¬ 
trum associated with the Gaussian part of the wavelet co¬ 
efficient distribution corresponds well to the Fourier power 
spectrum of |P| for I < 100 arcmin. The selection of the 
Gaussian part of the distribution is dependent on a param¬ 
eter in the scale-wise CVE algorithm which controls how 


restrictive the definition of an outlier is (Nguyen van yen 


et al. 20121. Nonetheless, the choice of this parameter is 


based on the symmetry of the separated “random” distribu¬ 
tions only. By comparing the total wavelet power spectrum 
of |VP| with the classical Fourier power spectrum of |P|, we 
see that the excess of power measured at intermediate scales 
might be associated with non-Gaussianities identified from 
the wavelet analysis. The amplitude of P and the ampli¬ 
tude of the gradient of P are two different tracers and they 
are not expected to produce the same power spectrum. The 
amplitude of the gradient of P traces changes in the polari¬ 
sation vector more strongly than the amplitude of P alone. 
Consequently, the excess of power measured at intermediate 
scales, partly displayed in Fig.[^ might be related to the ex¬ 
pected correlation between the magnetic field intensity and 
the electron density produced by compressive shocks. This 
measured excess of power suggests that the power spectrum 
of I VP| may trace fluctuations in the electron density as well 
as fluctuations caused by the Faraday rotation of polarised 
emission coming from a source localised behind the com¬ 
pressed magnetic field lines, whereas the power spectrum of 
|P| traces only fluctuations of the electron density. 

The power law index of the Gaussian power spectrum 
is shallower than that expected from a 3D Kolmogorov-like 
power spectrum (7 = 3.66) for a subsonic incompressible 
turbulent medium, however it is close to the index measured 
by Giardino et al. (2002 1 for the polarised intensity P at 2.4 


GHz from the Southern Galactic plane (7 = 2.37 ± 0.21). 
In the case of the non-Gaussian subset, many other physical 
processes can produce sharp changes in polarised data, such 
as outflows from massive stars and supernovae. We noted 
in section that features described as “double jumps” by 


Burkhart et al. (20121 are visible at different scales in the 


CGPS field. One of the two features framed by a green rect¬ 
angle (the upper right panel at Z = 45.7 arcmin) corresponds 
to the location of the supernova remnant SNR G84.2-0.8, 
which is clearly seen in the total intensity map. The “dou¬ 
ble jump” feature, which is itself part of the non-Gaussian 
subset, may be related to shocks produced by the supernova. 

On scales of 109 < I < 614 arcmin, we can see in Fig. 


12 (blue lines) that distributions become more symmetrical 


on larger scales. Since the distributions are normalised, it is 
reasonable to consider that the undersampling correction ap¬ 
plied to the 26-m Telescope data, which dominate on scales 
larger than ~ 80 arcmin, does not influence or bias distri¬ 
butions in a strong manner. The fitted power law with the 
Gaussian part of the distributions at lower scales could also 
be consistent with scales dominated by the DRAG 26-m, 
which could mean that non-Gaussian coefficients tend to ap¬ 
pear at smaller scales. According to Burkhart et al. (20121, 
the random distribution of Gaussian coefficients could be in¬ 
duced by subsonic turbulent fluctuations in electron density 
and magnetic field. On the other hand, non-Gaussian co¬ 
efficients may be associated with other physical processes, 
such as supersonic turbulence, gravitational collapses, mas¬ 
sive star outflows and supernovae. However, a clear separa¬ 
tion between lognormal distributions at smaller scales and 
Gaussian distributions at larger scales is hard to confirm. 
The lack of non-Gaussian contributions at larger scales could 
also be associated with the smaller number of statistically 
independent coefficients. A similar analysis realised on an 
extended range of scales could confirm if a real separation 
exists between distributions at smaller and larger scales. 


8 CONCLUSION 

We extend the calculation of |VP| to multiple scales us¬ 
ing a wavelet analysis formalism. The new technique shares 
similarities with the A-variance and the WTMM techniques 
used to characterise respectively the turbulence in molecular 
clouds and the multifractal nature of a surface or a medium. 
This approach can overcome the limitation of previous anal¬ 
yses which were only sensitive to the smallest scales. We 
show that fluctuations traced by IVPI exist at larger scales 
on data completed with lower spatial frequencies. Using the 
wavelet formalism, it is possible to measure the power spec¬ 
trum of |VP(Z,a;)| and evaluate the scaling behaviour of 
variations of the polarisation vector P in the Q-U plane. 
The scaling behaviour fallows a power law with 7 « 2.1. We 
measure a small drop in the spectrum between 80 < Z < 300 
arcmin. This drop corresponds well to the overlap in the 
u-v plane between the Effelsberg 100-m telescope and the 
DRAG 26-m telescope data. The undersampling presents in 
the 26-m telescope data has been identified as a source of 
unknown error in the data and could explain the measured 
drop of power. The wavelet analysis of | VP| also allows us to 
analyse the distribution of fluctuations in P as a function of 
angular scale. Distributions show higher skewness at smaller 
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scales than at larger scales. Separation of outliers contribut¬ 
ing to the tails of the distributions allows us to measure the 
power spectrum for the symmetrical part of the distribution. 
This power spectrum possesses a steeper power law with 
7 « 2.5. The spatial distribution of some outliers are part of 
correlated structures across angular scales, which trace the 
sharpest changes in the polarisation vector P in the field. 
Such higher intensity structures could be associated with 
compressive shocks of a supersonic turbulent medium. Fu¬ 
ture analysis applied over an extended range of scales and 
at higher Galactic latitude will provide a useful extension 
to the analysis presented here. Such analysis could confirm 
the appearance of a distinct type of fluctuation distribution 
at smaller scales as well as revealing the presence of high 
intensity structures at higher Galactic latitude or at higher 
angular scales. 
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